##########################################################################################################################################################
### Setup

### Source setup
source("code/1_modeling_setup.R")



##########################################################################################################################################################
### Analyze the effect of bicameralism itself by comparing Nebraska with similar states

### Figure out which states are a match for Nebraska in terms of the three nonbicameral SOP institutions in each year, 2010:2014

# Load institutions data, format and pare down so it's easy to look at
insts <- fread('data/stateleg-processvariables-final.csv',
               stringsAsFactors = FALSE, data.table = FALSE) %>%
  mutate(state_abb = state.abb[match(state, state.name)],
         state_year = paste(state_abb, year, sep='_'),
         chamber = str_extract(chamber, '^[A-Z]{1}')) %>%
  dplyr::select(state, year, chamber, veto_override_threshold, lineitemveto, lineitemveto_appropsonly) %>%
  dplyr::rename(lineitem=lineitemveto, veto_threshold=veto_override_threshold, lineitem_approps=lineitemveto_appropsonly) %>%
  arrange(state, chamber, year)

# Check SOP institutions values for Nebraska
filter(insts, state=='Nebraska') %>% # veto threshold is 0.6, has 1s for lineitemveto and lineitemveto_appropsonly
  .[1,] %>%
  dplyr::select(-year)

# Pull out states that match Nebraska's SOP institutions in each year
years <- 2010:2014
state_matches <- vector(mode='list', length=d(years))
names(state_matches) <- years
for(ii in 1:d(years)){
  df <- filter(insts, year==years[ii])
  ne <- filter(df, state=='Nebraska')
  res <- filter(df, veto_threshold==ne$veto_threshold & lineitem==ne$lineitem & state != 'Nebraska') %>%
    pull(state) %>%
    unique()
  state_matches[[ii]] <- res
}
state_matches <- state_matches %>% unlist() %>% unique() %>% .[!.=='North Carolina'] # NC is a miscode

# Make bicameral indicator
bill <- mutate(bill, bicameral = ifelse(state=='NE', 'Unicameral', 'Bicameral') %>% as.factor() %>% relevel(ref='Unicameral'))

# Filter out just Ohio (best SOP match and closest in terms of 1 party Republican control)
bicam <- filter(bill, state %in% c('NE', 'OH')) 

# Filter out OH, DE, IL, MD (all matches for veto_threshold and lineitem)
bicam2 <- filter(bill, state %in% c('NE', 'OH', 'DE', 'IL', 'MD'))


### Run models

# Reaching a vote
mods_unicam_vote <- vector(mode='list', length=2)
mods_unicam_vote[[1]] <- glm(voted_sles ~ bicameral + num_cs + as.factor(year), data=bicam)
mods_unicam_vote[[2]] <- glm(voted_sles ~ bicameral + num_cs + as.factor(year), data=bicam2)

# Becoming law
mods_unicam_law <- vector(mode='list', length=2)
mods_unicam_law[[1]] <- glm(law_sles ~ bicameral + num_cs + as.factor(year), data=bicam)
mods_unicam_law[[2]] <- glm(law_sles ~ bicameral + num_cs + as.factor(year), data=bicam2)


### Substantive effect sizes

# Vote prop
bicam_vote_sd <- bicam %>% pull(voted_sles) %>% sd(na.rm=TRUE)
bicam_vote_sd <- bicam2 %>% pull(voted_sles) %>% sd(na.rm=TRUE)

# Law prop
bicam_law_sd <- bicam %>% pull(law_sles) %>% sd(na.rm=TRUE)
bicam_law_sd <- bicam2 %>% pull(law_sles) %>% sd(na.rm=TRUE)

# Pr(vote)
bicam_effs_vote <- c(
  get_model_data(mods_unicam_vote[[1]], type='pred', terms='bicameral') %>% as.data.frame() %>% pull(predicted) %>% diff() %>% `/`(bicam_vote_sd),
  get_model_data(mods_unicam_vote[[2]], type='pred', terms='bicameral') %>% as.data.frame() %>% pull(predicted) %>% diff() %>% `/`(bicam_vote_sd)) %>%
  abs() %>%
  round(2)

# Pr(law)
bicam_effs_law <- c(
  get_model_data(mods_unicam_law[[1]], type='pred', terms='bicameral') %>% as.data.frame() %>% pull(predicted) %>% diff() %>% `/`(bicam_law_sd),
  get_model_data(mods_unicam_law[[2]], type='pred', terms='bicameral') %>% as.data.frame() %>% pull(predicted) %>% diff() %>% `/`(bicam_law_sd)) %>%
  abs() %>%
  round(2)


### Predicted effects (voted)

# Format estimates
coefs_unicam <- lapply(mods_unicam_vote, function(mod){
  dat <- plot_model(mod, type='pred', terms='bicameral')$data %>% as.data.frame()
  beta <- diff(dat$predicted)
  ci_lo <- beta - max(abs(dat$predicted - dat$conf.low))
  ci_hi <- beta + max(abs(dat$predicted - dat$conf.high))
  res <- data.frame(beta=beta, ci_lo=ci_lo, ci_hi=ci_hi)
  res
}) %>%
  do.call(rbind, .) %>%
  mutate(vars=c('OH', 'DE, IL, MD, OH'))

# Plot, format
plot_unicam_vote <- ggplot(coefs_unicam, aes(x=beta, y=vars, xmin=ci_lo, xmax=ci_hi)) +
  geom_pointrange(linewidth=2, size=1.5) +
  theme_bw() +
  labs(x='Effect of Bicameralism\nin Percentage Points', y='', title='Reaching a Vote') +
  geom_vline(xintercept=0, linetype=2, linewidth=2) +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        axis.text.x=element_text(angle=90, vjust=0.5),
        panel.border = element_blank(), 
        plot.title = element_text(hjust = 0.5)) +
  annotate(geom='text', x=0.06, y=2, label=paste(bicam_effs_vote[1], " SDs"), size=5) +
  annotate(geom='text', x=-0.045, y=1, label=paste(bicam_effs_vote[2], " SDs"), size=5)


### Predicted effects (law)

# Format estimates
coefs_unicam_law <- lapply(mods_unicam_law, function(mod){
  dat <- plot_model(mod, type='pred', terms='bicameral')$data %>% as.data.frame()
  beta <- diff(dat$predicted)
  ci_lo <- beta - max(abs(dat$predicted - dat$conf.low))
  ci_hi <- beta + max(abs(dat$predicted - dat$conf.high))
  res <- data.frame(beta=beta, ci_lo=ci_lo, ci_hi=ci_hi)
  res
}) %>%
  do.call(rbind, .) %>%
  mutate(vars=c('OH', 'DE, IL, MD, OH'))

# Plot, format
plot_unicam_law <- ggplot(coefs_unicam_law, aes(x=beta, y=vars, xmin=ci_lo, xmax=ci_hi)) +
  geom_pointrange(linewidth=2, size=1.5) +
  theme_bw() +
  labs(x='Effect of Bicameralism\nin Percentage Points', y='', title='Becoming Law') +
  geom_vline(xintercept=0, linetype=2, linewidth=2) +
  theme(text=element_text(size=16),
        axis.text=element_text(size=16),
        axis.text.x = element_text(angle=90, vjust=0.5),
        panel.border = element_blank(), 
        plot.title = element_text(hjust = 0.5)) +
  annotate(geom='text', x=-0.19, y=2, label=paste(bicam_effs_law[1], " SDs"), size=5) +
  annotate(geom='text', x=-0.095, y=1, label=paste(bicam_effs_law[2], " SDs"), size=5)


### Combine and save
plot_unicam <- plot_unicam_vote + plot_spacer() + plot_unicam_law + plot_layout(widths = c(4, 1 ,4))
ggsave("results/fig8_unicam_effects.pdf", plot=plot_unicam, width=15, height=6)
